sim.start.date = as.Date('2018-12-01')
# Keep a trailing 13 months. This is so we can compute 12 lagged differences.
price.path = copy(prices[date>=sim.start.date-months(1+lag.length), .(date, cpi, brent, wti, d.ln.wti, hh, d.ln.hh)])

# There's a gap between when FRED data ends (May) and futures data starts (July/Aug)
# Fill in those values
price.path = rbind(price.path, data.table(date=as.Date('2020-06-01')), fill=TRUE)
june.2020.cpi = cpi[DATE=='2020-05-01', CPIAUCSL] # not quite right but June CPI isn't available yet (until May 10th?)
june.2020.wti = wti[DATE=='2020-05-01', DCOILWTICO]*cpi.base/june.2020.cpi
june.2020.brent = brent[DATE=='2020-05-01', POILBREUSDM]*cpi.base/june.2020.cpi
june.2020.hh = hh[DATE=='2020-05-01', PNGASUSUSDM]*cpi.base/june.2020.cpi

price.path[date=='2020-06-01', cpi := june.2020.cpi]
price.path[date=='2020-06-01', hh := june.2020.hh]
price.path[date=='2020-06-01', wti := june.2020.wti]
price.path[date=='2020-06-01', brent := june.2020.brent]

price.path = rbind(price.path, data.table(date=as.Date('2020-07-01')), fill=TRUE)
price.path[date=='2020-07-01', wti := june.2020.wti]
price.path[date=='2020-07-01', brent := june.2020.brent]

# Add future dates
max.sim.date = as.Date('2050-12-01')
date.idx = seq.Date(price.path[,max(date)]+months(1), 
                    max.sim.date,
                    by='months')
price.path = rbind(price.path, data.table(date=date.idx), fill=TRUE)

# Add future price path
wti.fut = fread(data.dir%&%'/Energy Prices/Futures Strips/wti_futures.csv')
wti.fut[, date:=as.Date(date)]
wti.fut = wti.fut[,1:3]

brent.fut = fread(data.dir%&%'/Energy Prices/Futures Strips/brent_futures.csv')
brent.fut[, date:=as.Date(date)]
brent.fut = brent.fut[,1:3]

hh.fut = fread(data.dir%&%'/Energy Prices/Futures Strips/hh_futures.csv')
hh.fut[, date:=as.Date(date)]
hh.fut = hh.fut[,1:3]

price.path = merge(price.path, wti.fut[,.(date, WTI.futures.price)], by='date', all.x=TRUE)
price.path = merge(price.path, brent.fut[,.(date, brent_futures_price )], by='date', all.x=TRUE)
price.path = merge(price.path, hh.fut[,.(date, HH.futures.price)], by='date', all.x=TRUE)

price.path[, cpi := na.locf(cpi)] # replace NAs with last value

price.path[, WTI.futures.price := WTI.futures.price*cpi.base/cpi]
price.path[, brent_futures_price := brent_futures_price*cpi.base/cpi]
price.path[, HH.futures.price := HH.futures.price*cpi.base/cpi]

# Escalate last futures prices based on a linear regression (in first diffs)
# Project via a linear regression with MOY FEs
project.futures = function(varname, fut, date.range = price.path, log=FALSE) {
  dt.temp = copy(fut)
  if (log) dt.temp[, d.ln.p := c(NA, diff(log(dt.temp[[varname]])))]
  if (!log) dt.temp[, d.ln.p := c(NA, diff(dt.temp[[varname]]))]
  
  lm.fit = lm(d.ln.p ~ -1 + as.factor(month(date)), data=dt.temp)
  if (log) {
    pred.out = dt.temp[date==max(date)][[varname]]*
      exp(cumsum(predict(lm.fit, newdata=date.range[date>dt.temp[, max(date)]])))
  } else {
    pred.out = dt.temp[date==max(date)][[varname]]+
      cumsum(predict(lm.fit, newdata=date.range[date>dt.temp[, max(date)]]))
  }
  dt.temp = data.table(date=date.range[date>dt.temp[, max(date)], date], 
                       varname=pred.out)
  setnames(dt.temp, c('date',paste0(varname,'.projected')))
  return(dt.temp)
}
wti.project = project.futures(varname='WTI.futures.price', 
                              fut=price.path[!is.na(WTI.futures.price),.(date, WTI.futures.price)],
                              date.range=price.path, log=FALSE)
plot(WTI.futures.price ~ date , data=price.path, type='l', ylim=c(0,150))
lines(WTI.futures.price.projected ~ date, data=wti.project, col='red')

brent.project = project.futures(varname='brent_futures_price', 
                                fut=price.path[!is.na(brent_futures_price),.(date, brent_futures_price)],
                                date.range=price.path, log=FALSE)
plot(brent_futures_price ~ date , data=price.path, type='l', ylim=c(0,150))
lines(brent_futures_price.projected ~ date, data=brent.project, col='red')

hh.project = project.futures(varname='HH.futures.price', 
                             fut=price.path[!is.na(HH.futures.price),.(date, HH.futures.price)],
                             date.range=price.path, log=FALSE)
plot(HH.futures.price ~ date , data=price.path, type='l', ylim=c(0,5))
lines(HH.futures.price.projected ~ date, data=hh.project, col='red')

# The procedure yields prices that rise at average annual growth rates of 2.1 percent for WTI 
(wti.project[.N,2]/wti.project[1,2])^(1/
                                        (elapsed_months(wti.project[.N,date], wti.project[1,date])/12))
# and 2.5 percent for Henry Hub after 2030. 
(hh.project[.N,2]/hh.project[1,2])^(1/
                                      (elapsed_months(hh.project[.N,date], hh.project[1,date])/12))

price.path = merge(price.path, wti.project, by='date', all.x=TRUE)
price.path = merge(price.path, brent.project, by='date', all.x=TRUE)
price.path = merge(price.path, hh.project, by='date', all.x=TRUE)

# If we have the actual futures for a month, use that
price.path[!is.na(WTI.futures.price), wti := WTI.futures.price]
price.path[!is.na(brent_futures_price), brent := brent_futures_price]
price.path[!is.na(HH.futures.price), hh := HH.futures.price]

# If still missing, use projected futures
price.path[!is.na(WTI.futures.price.projected), wti := WTI.futures.price.projected]
price.path[!is.na(brent_futures_price.projected), brent := brent_futures_price.projected]
price.path[!is.na(HH.futures.price.projected), hh := HH.futures.price.projected]

# replace NAs with last value. only affects months without historical spot data OR futures.
price.path[, wti := na.locf(wti)] 
price.path[, brent := na.locf(brent)] 
price.path[, hh := na.locf(hh)] 

price.path = price.path[, .(date, wti, brent, hh)]
price.path[, d.ln.wti := c(NA, diff(log(wti)))]
price.path[, d.ln.hh := c(NA, diff(log(hh)))]
price.path[, d.ln.copper:=0] # need this if using the IV, for predict()
price.path[, spud_month:=date]
price.path[, moy:=month(date)]
